#include "synchronizability-calc.h"
#include "generator-algorithm.h"
#include <Eigen/Eigenvalues>
#include <Eigen/Dense>
#include <gtest/gtest.h>
#include <vector>
#include <algorithm>

namespace cccn
{

using namespace Eigen;
using namespace cccn;
using namespace boost;
using namespace std;

/**
 * \param adjm adjacency matrix
 * \return The lapican Matrix of A
 */ 
MatrixXd
GetLapicanMatrix (const MatrixXd& adjm)
  {
  MatrixXd diagm = adjm.rowwise ().sum ().asDiagonal ();
  return diagm - adjm;
  }

/**
 * \param g the graph object whose adjacency matrix is to be extracted as a MatrixXd object
 * \return adjacency matrix of g
 */
MatrixXd
GetAdjacencyMatrix (const Graph& g)
  {
  typedef graph_traits<BoostGraph_t>::edge_iterator edge_iterator;
  typedef graph_traits<BoostGraph_t>::vertex_descriptor vertex_descriptor;
  typedef graph_traits<BoostGraph_t>::edge_descriptor edge_descriptor;
  typedef property_map<BoostGraph_t, vertex_index_t>::type vertex_index_map_t;
  typedef property_traits<vertex_index_map_t>::value_type vertex_value_t;

  // get the numer of vertices in the graph
  const BoostGraph_t boost_g = g.GetBoostGraph ();
  uint n = num_vertices (boost_g);
  // create a n x n matrix, where n is the number of vertices in the graph
  MatrixXd adjm = MatrixXd::Zero (n, n);
  // loop through all the edges and set the matrix accordingly
  edge_iterator begin_e, end_e;
  tie (begin_e, end_e) = edges (boost_g);
  for (edge_iterator it = begin_e; it != end_e; it++)
    {
    // e is an edge in the graph
    edge_descriptor e = *it;
    // src and tgt are two vertices of the edge
    vertex_descriptor src = source (e, boost_g);
    vertex_descriptor tgt = target (e, boost_g);
    // src_i and tgt_i are the indices of the two vertices
    vertex_value_t src_i = get (vertex_index, boost_g, src);
    vertex_value_t tgt_i = get (vertex_index, boost_g, tgt);
    // set the adjacency matrix accordingly
    adjm (src_i, tgt_i) = true;
    adjm (tgt_i, src_i) = true;
    }

  return adjm;
  }

SynchronizabilityCalc::SynchronizabilityCalc (const Graph& g) : m_g (g)
  {
  MatrixXd adjm = GetAdjacencyMatrix (m_g);
  MatrixXd lapican = GetLapicanMatrix (adjm);

  SelfAdjointEigenSolver<MatrixXd> es (lapican);
  VectorXd ev_eigen = es.eigenvalues ();

  vector<double> ev;
  for (int i=0; i<ev_eigen.size (); i++)
    {
    ev.push_back (ev_eigen (i));
    }
  sort (ev.begin (), ev.end ());

  m_eigenvalue_2 = ev[1];
  m_eigenvalue_n = ev.back ();
  }

double 
SynchronizabilityCalc::GetSyncValue () const
  {
  return m_eigenvalue_2;
  }

double 
SynchronizabilityCalc::GetSyncRatio () const
  {
  return m_eigenvalue_2 / m_eigenvalue_n;
  }

} // namespace cccn

////////////////
//            //
// Test cases //
//            //
////////////////

using namespace cccn;

TEST (SynchronizabilityCalcUsageTest, ConstructorTest)
  {
  Graph g = GenerateTidyGraph_k3 (4);
  SynchronizabilityCalc sync (g);

  cout << "Lambda_2 = " << sync.GetSyncValue () << endl;
  cout << "Lambda_2 / Lambda_n = " << sync.GetSyncRatio () << endl;
  }

TEST (SynchronizabilityCalcInternalTest, GetAdjacencyMatrixFunctionTest)
  {
  Graph g = GenerateTidyGraph_k3 (4);
  MatrixXd adjm = GetAdjacencyMatrix (g);
  MatrixXd lapican = GetLapicanMatrix (adjm);

  SelfAdjointEigenSolver<MatrixXd> es (lapican);
  VectorXd eigenvalues = es.eigenvalues ();

  cout << "adjacency matrix = " << endl << adjm << endl;
  cout << "Lapican matrix = " << endl << lapican << endl;
  cout << "Eigenvalues = " << endl << eigenvalues << endl;
  }

TEST (SynchronizabilityCalcInternalTest, DISABLED_SelfAdjointEigenSolverTest)
  {
  MatrixXd X = MatrixXd::Random(5,5);
  MatrixXd A = X + X.transpose();
  cout << "Here is a random symmetric 5x5 matrix, A:" << endl << A << endl << endl;

  SelfAdjointEigenSolver<MatrixXd> es(A);
  cout << "The eigenvalues of A are:" << endl << es.eigenvalues() << endl;
  cout << "The matrix of eigenvectors, V, is:" << endl << es.eigenvectors() << endl << endl;

  double lambda = es.eigenvalues()[0];
  cout << "Consider the first eigenvalue, lambda = " << lambda << endl;
  VectorXd v = es.eigenvectors().col(0);
  cout << "If v is the corresponding eigenvector, then lambda * v = " << endl << lambda * v << endl;
  cout << "... and A * v = " << endl << A * v << endl << endl;

  MatrixXd D = es.eigenvalues().asDiagonal();
  MatrixXd V = es.eigenvectors();
  cout << "Finally, V * D * V^(-1) = " << endl << V * D * V.inverse() << endl;
  }
